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ABSTRACT 

High accuracy Monte Carlo simulation results for 1024*1024 Ising system 
with ferromagnetic impurity bonds are presented. Spin-spin correlation func- 
tion at a critical point is found to be numerically very close to that of a pure 
system. This is not trivial since a critical temperature for the system with 
impurities is almost two times lower than pure Ising T c . Small deviations 
from the pure behaviour contradict to some of the competing theories. 
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Influence of impurities on a critical behaviour has been a subject of numerous papers. 

For the simplest possible model - 2D Ising model this problem has been considered 
theoretically [1-10], experimentally [11, 12 1 and using computer simulations [13-18]. 

Most of the simulations were devoted to thermodynamic properties such as specific 
heat, magnetization and magnetic susceptibility. On the other hand, theories give direct 
predictions for spin-spin correlation function 

(S(0)S(r)} (1) 

where r is a distance between spins. 

We used cluster algorithm special purpose processor (SPP) J0|, [2(J to get accurate 
values of < S(0)S(r) > at a critical point. The SPP realizes in hardware Wolff cluster 
algorithm, and therefore does not suffer of critical slowing down. 

The SPP spends 375ns per one cluster spin. It also has a simple hardware for ex- 
tremely fast calculation of spin-spin correlation functions. Time necessary to get the 
correlation function for some r is equal to L 2 * 21ns, where L is a linear lattice size. 

We study the following model. Coupling constant J on each bond can take two 
positive values: J\ with probability p and Jq with probability 1 — p. For p = 0.5 duality 



relation |2T] shows that T c is equal to that of a pure model with all horizontal bonds 
equal to J\ and all vertical bonds equal to Jo- This greatly simplifies simulation data 
analysis, and so we used p = 0.5. 

Theoretical models use a small parameter 

g-piJo-J^ 2 (2) 

This parameter is connected with impurity induced length 

log/,-- (3) 
9 

To be able to notice deviations from the pure critical behaviour we should have 

k « L (4) 

To satisfy this condition we used L = 1024 and quite different values of J: Jq = 1 and 
Ji = 0.25. So in the simulations g is not very small. 

Theories deal with continuum limit and infinite lattice size. Simulations are conducted 
on a finite lattice with periodic boundary conditions. We calculated < S(0)S(r) > for 
spins, located along one lattice row. In this case distance r can take only integer values. 

Our simulations for the pure case [[H| showed that discrete lattice effects are significant 
for r < 8. Continuum theory can be applied for larger distances. But the finite size 
corrections for r > 8 are verv significant and should be taken into account explicitly. 



Pure Ising correlation function Co(r) for r/L — > has been calculated in [22 



, , 0.70338 

c °( r ) = -7174— ( 5 ) 

Continuum limit of (1) for the finite lattice with periodic boundary conditions c(r, L) has 
been obtained in [ f23| 1 

c ( r > L ) T74— ( 6 ) 



1/4 
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where 9 are Jacobi theta functions. For our purposes c(r,L) can be written in a simpler 
form 



1 + e"^ 4 |sin(f ) + cos(f )1 + e" 97r / 4 fcos(4r) - sin(^f 

C ^ L ) « 1 2 /- , L • ,o ,U (?) 

(sm(a) — e Z7r sm(3aj) 

where a = nr/L. The coefficient A(L) can be obtained using expression for co(r). 



Formula for c(r, L) is in exellent agreement with simulation results [2(| for the pure 
system. 

In Fig.l we show the ratio of computed < S(0)S(r) > to c(r, L) for L = 1024. 

To get mean values of < 5(0)5(r) > and standard deviations we used one thousand 
samples with different impurities distribution. For each sample all spins initially were 
pointing in the same direction. Two thousand Wolff clusters were flipped to thermalize 
spin distribution at critical temperature. Another eight thousand clusters were flipped 
to calculate mean values of < S(0)S(r) > for each sample. For L = 1024 one cluster flip 
at T c requires about O.lsec. 

Correlation functions for different r were measured for the same spin configurations. 

Error bars are determined mainly by different behaviour of < S(0)S(r) > for different 
impurity distributions and not by thermal fluctuations for a given sample. 

Deviations of < S(0)S(r) > from c(r, L) at r < 8 are due to the discrete lattice effects. 
In fact, at these distances correlation function of the impure system is extremely close to 
the pure correlation function. For r = 1 difference between them is of the order of 10 -3 . 
At first sight it seems to be quite natural, because r << U and < S(0)5(0) >= 1. But 
the critical temperatures for pure and impure cases are different, and continuum theory 
cannot exclude strong renormalization of < 5(0) 5(1) >. 

On the other hand, for r > 8 pure correlation function practically coincides with 
c(r, L). So, Fig.l shows that impurities decrease spin-spin correlations at large r. Again, 
this is not trivial because of the difference of critical temperatures for pure and impure 
cases. 

This result contradicts the renormalization group calculation of DD |l] for the spin- 
snin correlation function, which would instead show increase if spin correlations. - see 



also the discussion in the first paper of [ 13 1 where the numerical simulation problem has 
been defined. 



The calculations of DD were based on the interacting fermion model |2J, £5(] repre- 
senting the 2D Ising model with impurities p6| . The replica componets trick has been 
used, for the fermion fields renormalization group calculation, which gave a new critical 
behavior for the correlation length and of the specific heat, as compared to the Ising 
model on a perfect lattice. The result of [J25| has later been confirmed by calculations 
using other approaches, see [2-7], as well as by the computer simulations [13]. (See how- 
ever [8,9] which claim a finite specific heat at T c ). In these calculations, of the energy 
operator related quantities, the replica components number N, to be taken eventually 
equal to 0, plays a minor role, that of the diagram counting, both in the original calcula- 
tion in [26] using fermions, as well as in other approaches (bosonization, spin-component 
replicas) [2-7] . 

On the other hand, in the calculation of the spin-spin correlation function in [27] 
(see also [1]) the fermion replica components did play a dynamical role, while the other 
approaches, dealing more directly with the spin operators, employed replicas for diagram 
counting only. 

Our present numerical result for the spin-spin function shows that fully dynamical 
effects due to replicas of the fermion fields lead to a wrong result, for this particular 
problem. The precise reason for this remains unclear. 

The approaches based on bosonization gave the spin-spin function to be exactly same 
(in the scaling limit) as that of the pure Ising model. This is also not supported by our 
numerical result. 

The diviations of (SS) in Fig.l should be compared then with the standard renor- 
malization group applied directly to the spin operators. The result of [2,5,6] was that the 
coefficient of the first order renormalization group equation vanishes, in the iV — > limit 
(N now is the spin operator replicas). No one had actually calculated the second order 
coefficient (let us denote it C2), but assuming it is non-vanishing the result is standard. 
As the renormalization group equation one has: 

= -^ 2(x) (8) 

Here 

(g(O)S(r)) 

HX} (S(0)S(r)) o W 

(S(0)S(r)) ~ co{r) ~ (10) 

x = log(r) (11) 



The renormalization group coupling constant g(x) is [26|: 



7V 



Then from (|S|) one finds: 



and 

^(0)%)>-^I^P(-T^) ( 14 ) 

7T 

So it is not exactly the pure Ising model behavior as there is a cross-over in the 
amplitude: 

T(x) ~ 1, x <C Xi 

« exp( — po), x>Xj (15) 

Xj = vr/4(7o = log(Zj) is the impurities driven cross-over scale. 

In the cross-over region (SS) / (SS)o = T(x) is decreasing. This is what we observe 
in Fig.l. 

Suppose that the finite lattice correlation function renormalization is given by (13) 
as well. Then we can estimate C 2 . Value of go, found from the magnetic sucseptibility, 
magnetization and specific heat data [13,28], is somewhere near 0.3. According to Fig.l 
C*2 turns out to be quite small, of the order of 1CT 2 . This also agrees with the magnetic 
susceptibility and magnetization data [28]. 

More specifically, the estimate by the slope of the plot in Fig.l for r relatively small, 
8 < r < 32, where the finite size corrections should still be insignificant, gives C 2 ~ 0.03. 
Clearly, this is an order of magnitude estimate. (Here, by finite size corrections, we mean 
in fact the small corrections to them due to the extra factor T(x), eqs.(13),(14), in the 
correlation function of the impure model, as compared to the pure one). 



We remark finally that previous measurements of the magnetization in [13] has not 
shown difference with the pure model because the precision was not sufficient. In fact, 
the variation we observe for (SS) / (SS) corresponds, for small 4g /nx, to an effective 
change in the magnetization index on the value Cig\j2, which is from 0.125 to ~ 0.126. 



This is just beyond the accuracy in |L3 . 
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Figure Captions 



Fig.l. Ratio of < S(0)S(r) to c(r,L). Solid line connects data points for the system 
with impurities at (1/T C ) = .8070519. Dashed line shows data for the pure Ising model 
at (1/T C ) = .4406868. 



